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ABSTRACT 

The inner regions of barred galaxies contain substructures such as off-axis shocks, nuclear rings, and nuclear 
spirals. These substructure may affect star formation, and control the activity of a central black hole (BH) 
by determining the mass inflow rate. We investigate the formation and properties of such substructures using 
high-resolution, grid-based hydrodynamic simulations. The gaseous medium is assumed to be infinitesimally- 
thin, isothermal, and non-self-gravitating. The stars and dark matter are represented by a static gravitational 
potential with four components: a stellar disk, the bulge, a central BH, and the bar. To investigate various 
galactic environments, we vary the gas sound speed c s as well as the mass of the central BH Mbh. Once the 
flow has reached a quasi-steady state, off-axis shocks tend to move closer to the bar major axis as c s increases. 
Nuclear rings shrink in size with increasing c s , but are independent of Mbh, suggesting that ring position is not 
determined by the Lindblad resonances. Rings in low-c^ models are narrow since they are occupied largely by 
gas on X2 -orbits and well decoupled from nuclear spirals, while they become broad because of large thermal 
perturbations in high-c^ models. Nuclear spirals persist only when either c s is small or M B h is large; they would 
otherwise be destroyed completely by the ring material on eccentric orbits. The shape and strength of nuclear 
spirals depend sensitively on c s and M B h such that they are leading if both c s and M B h are small, weak trailing if 
c s is small and Mbh is large, and strong trailing if both c s and Mbh are large. While the mass inflow rate toward 
the nucleus is quite small in low-c^ models because of the presence of a narrow nuclear ring, it becomes larger 
than 0.01 M yr -1 when c s is large, providing a potential explanation of nuclear activity in Seyfert galaxies. 
Subject headings: hydrodynamics — galaxies: ISM — galaxies: kinematics and dynamics — galaxies: nuclei 
— galaxies: spiral — ISM: general — shock waves 



1. INTRODUCTION 

Stellar bars play an important role in the dynamical 
evolution of gas in galaxies. By introducing a non- 
axisymmetric torque, they produce interesting morphological 
substructures in the gaseous medium, including a pair of 
dust lanes at the leading side of the bar, a nuclear ring 
near the center, and nuc le ar spirals insi d e the r ing (e.g., 
Sanders & Huntlevl [T97g; |Eoberts_eIiD U979|; ISchwarzl 
T98T1: Iva n Albada & Roberts 1981; Athanassou lj U992b|; 
Piner etal.1 119951: iButa&Combel 119961: iMartini et all 
2003a b; Martinez- Valp uesta et al.ll2006|) . They also tra nsport 
gas inward which can trigge r starb ursts in the rings (e.g ., lButal 
19861: iGar cia-Barr eto et ap [199JJ; [gefleL & Shlosman 19941: 
Barth et al . 1995; Maoz et al.ll200Tl: iMazzuca et al.ll2008l) and 
if the mass inflow extends all the way to the c enter, they may 
help p o wer active galactic nu cl ei (AG N ) (e.g.,[S hlosman ^et al] 
1990; Regan & Mulc haevl fl999l: iKnapen et al.1 l2000l: 
Laurikainen et al. 2004; van de Ven & Fathi 2010). 

Since bar substructures represent a nonlinear response 
of the gas to a non-axisymmetric gravitational potential, 
their formation and evolution is best studied using direct 
numerical simulations Q There have been a number of 
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1 Englmaier & Shlosman ( 2000) argued that physical properties of nu- 
clear spirals can be explained by the linear density-wave theories (see also 
Macieiewski 2004a). 



numerical studies on the gas dynamics in barred galax- 
ies. Based on the numerical scheme employed, they can 
be categorized largely into two groups: (1) those using 
a smoothed particle hydrodynamics (SPH) techniqu e (e.g. , 
lEnglmaier & Gerhard! 19971: IFa tsis & Atha nassoulal [2000; 
lAnn & Thakuri 120051: iThakur et al.1 2009) and (2) those 
using a gr i d-base d algorithm (e.g., lAt hanassoula 1992b; 
Pineretal. 1995; Macieiewski et al. 2002; Macieiewski 
I2004fi iRegan &Teubenll2003L 120041) . The numerical results 
from these two approaches do not always agree with each 
other, at least quantitatively, even i f the mod e l para meters 
are almost identical. For instance, iPiner et all (119951) using 
the CMHOG code on a cylindrical grid reported that the gas 
near the corotation regions exhibits complex density features 
resulting from Rayleigh-Taylor and/or Kelvin-Helmholtz in- 
stabilities, while these structures are absent in the SPH sim- 
ulations. In addition, overall shapes and structures of dust 
lanes and nuclear rings from CMHOG simulations are differ- 
ent from SPH results. 

Some differences in the numerical results may be at- 
tributable to relatively large numerical diffusion of a standard 
SPH method an d its inability to handle sharp discontinuities 
accur ately (e.g., lAgertz et al.1 120071 : iPricd 120081: iRead et al.1 
2010). However, after adopting and thoroughly testing the 
CMHOG code as part of this work, we have found it con- 
tained a serious bug in the way the gravitational forces due 
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to the bar are added to the hydrodynamical equations. Thus, 
some of the discrepancies in the flows computed by CMHOG 
and other codes are likely due to this bu g. We discuss this b ug 
and its affect on the results reported in lPiner et al.l (Il995|) in 
Section 2.2. 

In this paper, we revisit the gas dynamics in barred galaxies 
using a corrected version of the CMHOG code. Our objec- 
tives are three-fo ld. First, we wish to remedy the errors in 
iPiner et all (Il995|) . and to compute the formation of bar sub- 
structures with an accurate shock-capturing grid code with 
the correct bar potential. Second, the morphology, shape, 
and strength of the bar substructures are likely to depend on 
the gas sound speed a nd the shape of the underlyin g gravi- 
tational potential (e.g., Englmaier & Shlosman 2000). Thus, 
we report new models in which we include a central black 
hole (BH) that greatly affects the gravitational potential in 
the central regions, and we vary both the BH mass as well 
as the sound speed to explore the dynamics in various galac- 
tic conditions. Third, we exploit advances in computational 
resources to compute models that have more tha n an order 
of ma gnitude higher resolution than the models in lPiner et al.l 
(1995), with a grid resolution of 0. 13 pc in the central regions. 
This allows us to resolve details in the flow in the nuclear re- 
gions, in particular the formation of nuclear rings and nuclear 
spirals. 

According to the most widely accepted theory, a nuclear 
ring forms near the inner Lindblad resonance (ILR) when 
there is only one ILR, as the gas outside (inside) ILR loses 
(gains) angular momentum and accumulates there, while it 
forms in between th e inner ILR and outer I LR when there 
are two ILRs (e.g., j S Mosman et alJ 1 19901: ICombesI 119961: 
lButa&Combeslll996l) . On the other hand. iRegan & Teubenl 
d2003h argued that the ring formation is more deeply related 
to the existence of x^ -orbits rather than the ILRs. But, the ar- 
guments relying on either ILR or x^ -orbits do not take into ac- 
count the effect of thermal pressure. Therefore, it is important 
to explore to what extent the concepts of ILR or x^ -orbits are 
valid in describing nuclear rings, especially when the sound 
speed is large. 

The formation, shape, and nature of nuclear spirals that may 
channel the gas to the galaxy centers are also not well under- 
stood. Observations using the Hubble Space Telescope indi- 
cate that galaxies having nucle ar dust spirals are quite com- 
mon (e.g., Ma rtini et al. 2003a : b). While most of such spirals 
are trailing, a few galaxies including NGC 1241 and NGC 
6902 reportedly po ssess leading nuclear spirals ( Dia z et al.l 
120031: iGro sbol 2003). Although the linear theory suggests that 
leading spirals are expected when there are two ILRs (e.g., 
iMaci eiewski 2004a), they are absent in the numerical models 
of lPiner et al.1 (Il995|) computed with the C MHOG code, while 
the SPH models of lAnn & ThakuH (120051) with self-gravity do 
form leading spirals. The SPH models suffer from poor spa- 
tial resolution at the nuclear regions as most particles gather 
around the rings. By running high-resolution simulations with 
a corrected version of CMHOG, we can clarify the issues of 
the nuclear spiral formation and related mass inflow rates to 
the galaxy center. 

We in this work treat gaseous disks as being two- 
dimensional, isothermal, non-self-gravitating, and unmagne- 
tized, which introduces a few caveats that need be noted from 
the outset. By considering an infinitesimally-thin disk, we ig- 
nore gas motions and associated dynamics along the direction 
perpendicular to the disk. By imposing a point symmetry rela- 
tive to the galaxy center, our models do not allow for the exis- 
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a The BH mass is varied with time as Mbh(0 = ^bh(0) + 
J^Mit^dt' assuming that all the inflowing mass is added 
to the central BH. 

tence of odd-m modes, although this appears reasonable since 
m = 2 modes dominate in the problems involving a galactic 
bar. In addition, we are unable to capture the potential conse- 
quences of gaseous self-gravity and magnetic stress that may 
not only cause fragmentation of high-density nuclear rings but 
also affect mass inflow rates to the galaxy center. Neverthe- 
less, these idealized models are useful to isolate the effects 
of the gas sound speed and the mass of a central BH on the 
formation of bar substructures and mass inflows. Also, these 
models allow us to correct the results of previous CMHOG 
calculations with incorrect bar forces. 

This paper is organized as follows. In Section 2, we de- 
scribe the galaxy model, model parameters, and our numerical 
methods. In Section 3, we present the results of simulations 
for off-axis shocks and nuclear rings. The detailed properties 
of nuclear spirals are presented in Section 4. In Section 5, we 
study the mass inflow rates through the inner boundary ob- 
tained from our simulations. In Section 6, we conclude with a 
summary and discussion of our results and their astronomical 
implications. 

2. MODELS AND METHODS 

We consider a uniform, isothermal, infinitesimally-thin, and 
non-self-gravitating gas disk orbiting in a gravitational poten- 
tial $ext arising from various components of a barred galaxy. 
The bar is assumed to rotate about the galaxy center with a 
fixed pattern speed = fi^z. Therefore, it is advantageous to 
solve the dynamical equations in cylindrical polar coordinates 
(R, (j)) corotating with the bar in the z = plane. The equations 
of ideal hydrodynamics in this rotating frame are 

^+u-vjs = -EV-u, (1) 

f_ + u .vju = -c?-- V$ext + - 2Q h X u, (2) 

where E, u, and c s denote the surface density, velocity in the 
rotating frame, and the sound speed in the gas, respectively. 
The third and fourth terms in the right hand side of equation 
© represent the Coriolis and centrifugal forces, respectively, 
arising from the coordinate transformation from the inertial 
to rotating frames. The velocity v in the inertial frame is ob- 
tained from v = u + Rflb4>. In order to focus on the bar-driven 
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gas dynamics, we do not consider star formation and the as- 
sociated gas recycling in the present work. 

The external gravitational potential <£ ext consists of four 
components: an axisymmetric stellar disk, a spherical bulge, 
a non-axisymmetric bar, and a central supermassive BH. Ap- 
pendix El describes the specific potential model we employ 
for each component of the galaxy. The bar pattern speed is 
taken to be = 33 km s -1 kpc -1 . Without a centra l BH, o ur 
galax y model i s simil ar to those in Athanass oulal (Il992albh 
and iPiner et al.l (Il995|) . The presence of a BH allows us to 
explore the effect of central mass concentration on the forma- 
tion o f nuclear spirals (e.g., Macieiewski 2004b HThakur et all 
l2QQ9h . 

2.1. Models 

The real interstellar gas is multiphase and turbulent, with 
temperatures differing by a few orders of magnitude (e.g. , 
Field . Goldsmith. & Habingl [T969I: iMcKee & Ostriker] fl977l 
2007). For simplicity, we model this highly inhomogeneous 
gas using an isothermal equation of state with an effective 
sound speed c s that includes a contribution due to turbulent 
motions. We have calculated 15 different models in which we 
vary both c s and the initial mass Mbh(0) of the central BH as 
parameters. Table [T] lists the properties of each calculation. 
The sound speed is chosen to vary between 5 and 20 km s -1 . 
Models with a postfix bhO do n ot initially possess a central 
BH, and are similar to those in IPiner et al.l (Q995). Models 
with the postfix bh7 or bh8 have a BH with mass M B h(0) = 
4 x 10 7 M and 4 x 10 8 M^ resp ectively; they are analogou s 
to models in lMacieiewskil (I2004bl) and lAnn & Thakuri {2005). 
For most models, we fix the BH mass to its initial value, but 
we also consider three additional models (cs05bh0t, cs20bh0t, 
and cs20bh7t) in which the BH mass is varied with time ac- 
cording to M B h(0 = M B h(0)+ ^M(t')dt where M is the mass 
inflow rate across the inner boundary (see below). These time- 
varying Mbh models allow us to study the effect of BH growth 
due to the gas accretion on bar substructures. 

Figure Q] plots the net circular rotation curves together with 
a contribution from each component when M B h = 4 x 10 8 M . 
The solid and dashed lines are along the bar major and minor 
axes, respectively. The circular velocity is almost flat at ~ 
200 km s -1 in the outer parts. Without the BH, the rotation 
curve v c would rise linearly with R close to the center, but 
the presence of the BH results in v c oc R~ 1 ^ 2 for R < 0.1 kpc. 
This rapid increase of v c will render the gaseous orbits in the 
very central regions highly resistant to pressure perturbations, 
resulting in smaller mass inflow rates than the cases without 
it, as we show below. 

Figure [2] shows the characteristic angular frequencies, 
Q — k/2, tt, and Q + k/2 along the bar major and minor 
axes as solid and dashed lines, respectivelyQ Here, Q 2 = 
R~ l d$ QXt /dR and k 1 = R~ 3 d(R 4 Q 2 )/dR denote the angular 
and epicyclic frequencies, respectively. The horizontal dot- 
ted line in each panel represents the bar patten speed of Vt^ = 
33 km s -1 kpc -1 , with the corotation resonance (CR) located 
at Rqr = 6 kpc for all models. For bhO models with no BH, the 
Q-k/2 curve peaks at R max = 0.53 kpc and is equal to fl^ at 

2 In the presence of a non-axisymmetric bar potential, the concepts of Q as 
an angular frequency and k as a radial frequency do not apply strictly since 
closed orbits are in general non-circular. In our models, however, the bar 
potential is nearly axisymmetric at R < 1 kpc, so that Q and k measure the 
actual frequencies reasonably well in the central parts. 
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Fig. 1 . — Rotational velocity of each component of the model galaxy with 
a central BH of Mbh = 4 x 10 8 Mq . The solid and dashed lines are for along 
the bar major and minor axes, respectively. Note that the effect of the BH 
is almost negligible at R > 1 kpc, while it dominates the total gravitational 
potential at R ^ 0. 1 kpc. 



the two ILR with radii of 7?hlr = 0. 19 kpc and 7? ilr ~ 2 kpc. 
Because the rotation curve rises steeply toward the center, bh7 
and bh8 models have only a single ILR at 7?ilr ~ 2 kpc. The 
Q-k/2 curve in bh7 models attains its local maximum and 
minimum at R max = 0.53 kpc and R m [ n = 0.19 kpc, respectively. 
In bh8 models, on the other hand, it increases monotonically 
with decreasing R since the BH dominates the gravitational 
potential. We will show in Section [4] that the shape of nuclear 
spirals depends critica lly on the sign of d{Vt- n/2)/dR (e.g., 
lButa&Combell996l) . 

2.2. Numerical Methods 

To solve equations (Q]) and 0, we use the two-d imensional 
grid-b ased code CMHOG in cylindrical geometry (IPiner et alJ 
|1995|) . CMHOG implements the piece wise parabolic method 
in its Lagrangian remap formulation (Colell a~& Woodward! 
119841) . which is third-order accurate in space and has very lit- 
tle numerical diffusion (viscosity). All the runs are carried out 
in a frame corotating with a bar whose major axis is aligned 
along the y-axis (i.e., <p = ±ir/2), so that the bar potential re- 
mains stationary in the simulation domain. By assuming a re- 
flection symmetry with respect to the galaxy center, the simu- 
lations were performed on a half-plane with -tt/2 < (j) < tt/2 
constructed by making a cut along the bar major axis. 

As me ntioned in S e ction [T] the original version of CMHOG 
used by IPiner et all (Il995|) contained a serious bug in the 
way the gravitational forces were added to the hydrodynamic 
equations. The CMHOG code places a bar potential 3>bar 
with the major axis aligned along the x-axis, calculate the 
bar forces f x = -dQ^/dx and f y = -d^^/dy, and then 
transforms them into cylindrical coordinates. In the orig- 
inal version of the code, the incorrect transform relations 
f R = f x cos (j) + f y sin (j) and f^ = f x sin (j)-f y cos (j) were used. 
In fact, the correct transformation rule for the azimuthal force 
should be f^ = -f x sin 4> + f y cos <j). With the sign of the az- 
imuthal force reversed (but the radial force correct), the flows 
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Fig. 2. — Angular frequencies of galaxy models with different BH masses. 
The solid and dashed lines represent Q — k/2 (rightmost curves), Q (middle 
curves), and Q + k/2 (leftmost curves) along the bar major and minor axes, 
respectively. The dotted lines denote the bar pattern speed Q^. (a) Models 
without BH have two ILRs at Rulr = 0.19 kpc and /?oilr ~ 2 kpc, with the 
maximum of the Q-k/2 curve occurring at 7? max = 0.53 kpc. (b) Models 
with M B h = 4 x 10 7 M have a single ILR at Ri LR « 2 kpc, with the local 
maximum and minimum of the Q-k/2 curve occurring at R m ax = 0.53 kpc 
and R min = 0. 19 kpc. (c) Models with M B h = 4 x 10 8 M have a single ILR at 
Rjlr ~ 2 kpc, with d(Q-K/2)/dR < in the nuclear regions with R < 1 kpc. 

in models computed using the original CMHOG code behave 
as if the bar potential were aligned parallel to the v-axis, but 
with forces quite different from the intended ones. Other than 
these force transformations, the complex hydrodynamic algo- 
rithms in CMHOG are all imple mented correctly, a nd were 
well tested in the original paper by lPiner et al.l (Il995|) . There- 
fore, previous numerical studies based on CMHOG should 
remain valid as long as they did not adopt the incorrect trans - 
formations of the bar forces inherited fromlPiner et al.l (119951) . 
Unfortunately, the results of Pine r et all (Il995|) were compro- 
mised by a trivial sign error in the coordinate transform rela- 
tions for the bar forces. 

Figure 3 compares the results for a typical simulation run 
with the original and corrected version of CMHO G. For this 
test, the grid resolution is taken identical to that in Pine r et all 



(119951 ) with 251 and 154 zones in the radial and azimuthal 
directions, respectively. The figure plots the logarithm of the 
surface density at t = 300 Myr. Several differences are ap- 
parent. For example, the gas around the corotation resonance 
(CR) at R = 6 kpc in the left panel is largely evacuated and 
has corrugated streams linked to the ends of the bar major 
axis, whereas the CR region is relatively featureless in the 
right panel. The nuclear ring in the left panel is fairly smooth, 
while it is quite clumpy in the right panel. Most importantly, 
the very central region inside the ring is almost unperturbed 
in the left panel, while the right panel shows spiral structures 
in the central region. These differences suggest that the orig- 
inal CMHOG code is unable to properly model the flow in 
the central regions, especially weak nuclear spirals. We have 
also run the same model using other gri d-based codes a dopt- 
ing Cartesia n coordinates, such as TVD (I Kim et al.lll999|) and 
ANTARES (lYuan & Yenll2005f) as w ell as the particle-based 
GADGET code (Spring el et alJl2001h . in all of which the bar 
forces are calculated by taking finite differences of <l>bar rather 
than using f x and f y directly. The new version of the CMHOG 
code used in this work gives results which are much more sim- 
ilar to the results of these other codes, which gives us further 
confidence that the gravitational forces due to the bar are now 
being treated correctly. 

To resolve the central regions accurately, we set up a non- 
uniform, logarithmically- spaced cylindrical grid with 1024 
radial zones extending from 0.02 kpc to 16 kpc and 480 az- 
imuthal zones covering the half-plane. This makes the zones 
approximately square- shaped throughout the grid (i.e., AR = 
RAcj)). The resulting grid spacing is AR = 0.13, 6, and 100 pc 
at the inner radial boundary, R = 1 kpc where nuclear rings 
typically form, and the outer radial boundary, respectively. 
This increases the resolution in the inner regions by over an 
or der of magnitude, in comparison to the models presented 
in iPiner et aD (Il995|) . This level of grid resolution is cru- 
cial to resolve nuclear spiral structures within R = 1 kpc. We 
use outflow and continuous boundary conditions at the inner 
and outer radial boundaries, respectively, while the azimuthal 
boundaries are periodic. The gas crossing the inner boundary 
is considered lost from the simulation domain. We keep track 
of the total mass crossing the inner boundary in order to study 
the mass inflow rates into the galactic nucleus. 

Each model starts from a uniform disk with surface density 
Ho = 10 M pc~ 2 that is rotating in force balance with an ax- 
isymmetric gravitational potential without a bar. In order to 
avoid strong transients in the fluid flow caused by a sudden 
introduction of the bar, we slowly introduce the bar potential 
over one bar revolution time of 27r/^ = 186 Myr. This is ac- 
complished by increasing the bar central density pbar linearly 
with time and decreasing the bulge central density pt>ub while 
keeping the net central density pbar + Pbui fixed. This ensures 
that the shape of the total gravitational potential <£ e xt, when 
averaged along the azimuthal direction, is unchanged with 
time. All the models are run until 500 Myr. This corresponds 
to 1.2 x 10 4 and 10 orbits at the inner and outer radial bound- 
aries, respectively, for bh8 models with M B h = 4 x 10 8 M©. 

3. RESULTS 

We take Model cs05bh0 with c s = 5 km s -1 and no BH as 
our standard model. The overall evolution of other models 
with different c s and Mbh are qualitatively similar, although 
the properties of the nuclear features that form differ consid- 
erably from model to model. In this section, we first describe 
the evolution of our standard model, and then present the dif- 
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Fig. 3. — Logarithm of the gas surface density at t = 300 Myr from a test run using (a) the original CMHOG code and (b) the corrected version used in this 
work. A cylindrical grid with 251 x 154 zones is used. Compared to the left panel, gas in the right panel is relatively featureless in the corotation region at 
R = 6 kpc, has a more clumpy ring, and harbors nuclear spirals in the central region. 



ferences in the off-axis shocks and nuclear rings caused by 
differing c s and Mbh- The properties of nuclear spirals will be 
given in Section |4l 

3.1. Overall Evolution 

Figure [4] plots snapshots of the logarithm of the gas den- 
sity at a few selected epochs in the inner regions of Model 
cs05bh0. The bar is oriented vertically along the y-axis, and 
the gas is rotating in the counterclockwise direction relative 
to the bar. The images extend to 6 kpc on either side of the 
center, corresponding to the CR radi us, outside of wh ich the 
gas remains almost unperturbec0. iPiner et all (Il995|) found 
that the CR regions exhibit time-dependent flow structures, as 
reproduced in Figure [3^. On the other hand, Figure 0] shows 
that the CR regions in our simulations are quite stable and 
exhibit only at late time low-amplitude wavelike features en- 
trained by the dense gas located at the bar ends, similar to the 
results of SPH simulations ( e.g., Englmaier & Gerhard 1997; 
iPatsis & Athanasso ula 2000 ). This indica tes that the compli- 
cated structures seen in Piner et al. (1995) were likely an arti- 
fact of the errors in their force transformation. 

A striking feature of each of the snapshots shown in Figure 
Slat times greater than 120 Myr is large amplitude oscillations 
in the density in the ring and dust lanes. These fe atures are 
also seen in the simulations of Wada & Koda ( 2004) for spiral 
shocks, and are attributed to the "wiggle instability". As we 
will discuss below, this shock instability appears to be caused 
by vorticity generation in curved shocks. 

An introduction of the non-axisymmetric bar potential in- 
duces perturbations on the gas orbits, causing them to devi- 
ate from circular trajectories. The gas density increases (de- 
creases) in regions where neighboring orbits come close to- 
gether (diverge). When t = 60 Myr, the overdense regions 
are preferentially located downstream from the bar major axis 
(Fig. |4k). At this time, the overdensity produced by orbit 

3 The non-axisymmetric bar potential we adopt is very weak at R > Rqr. 



crowding is largest at R ~ 3 kpc. Since the perturbing force 
by the bar is weak inside R~ 1 kpc, the overdensity there is 
correspondingly small and not readily discernible. Over time, 
the overdense regions become narrower and sharper as the bar 
amplitude grows and eventually develop into shock fronts at 
around t = 100 Myr. In what follows, we term these narrow 
shocks off-axis shocks. 

A nuclear ring is beginning to shape at this time, as well. 
To illustrate the formation of nuclear rings in our models, we 
plot as solid lines in Figure \5\ instantaneous streamlines of the 
gas that starts from Point A marked at (x,y) = (1.5,-2.7) kpc 
in Model cs05bh0 and from (x,y) = (1.4,-2.5) kpc in Model 
cs20bh0 with c s = 20 km s" 1 and no BH at t = 100 Myr. The 
two dotted circles in Figure [5^ indicate the inner and outer 
ILRs at Rjjl R = 0.19 kpc and R 0JLR = 2 kpc. Note that the 
thick lines representing the overdense ridges in both models 
directly cross the outer ILR. The changes of the azimuthal 
and radial velocities along the streamlines are shown in Figure 
[5j),c, where the dotted line indicates the equilibrium rotation 
curve of the model galaxy with no BH. On emerging from 
the overdense region (Point A), the gas moves radially inward 
on its epicycle orbit and increases (decreases) its azimuthal 
(radial) velocity due to the Coriolis force. It reaches Point 
B closest to the center when it attains vr = and largest v^. 
After this point, it moves radially outward, decreasing until 
it hits the off-axis shock at Point C. The gas loses a significant 
amount of angular momentum at the shock and begins to fall 
in. In addition, the shocked gas is swept by other shocked gas 
flowing from the bar end regions along the shocks. Note that 
the shape of the off-axis shocks shown in Figured coincides 
with the gas streamline from Points C to D, indicating that all 
the gas after crossing the shocks moves radially in along the 
shock fronts in the developing stage of the nuclear rings. 

As the shocked gas moves along the shock fronts from Point 
C, it gradually rotates faster again. When the azimuthal veloc- 
ity of the gas is increased to the level comparable to the equi- 
librium circular velocity at some radius R, it begins to follow 



6 



Kim et al. 




x (kpc) x (kpc) 



Fig. 4. — Snapshots of the logarithms of the gas surface density of Model cs05bh0. The bar is oriented vertically along the v-axis and remains stationary. The 
gas inside the CR is rotating in the counterclockwise direction. In (f), the dotted curves aligned vertically represent the x\ -orbits that cut the x-axis at x c = 0.4, 
0.8, 1.2, 1.6 kpc from inside to outside, while the solid curves aligned horizontally plot the X2 -orbits with x c = 0.2, 0.6, 1.0, 1.4 kpc. Clumpy structures in (c)-(e) 
are produced by vortex generation at the curved shocks. See text and Figure [6] for detail. 
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Fig. 5. — (a) Instantaneous streamlines of gas that starts from Point 
A (x,y) = (1.5,-2.7) kpc in Model cs05bh0 and (1.4,-2.5) kpc in Model 
cs20bh0 at t = 100 Myr. The thick green and orange curves represent the 
overdense ridges in Models cs05bh0 and cs20bh0, respectively. The two dot- 
ted circles indicates Rulr = 0.19 kpc and /?oilr = 2 kpc. In Model cs05bh0, 
the gas reaches Point B closest to the galaxy center, is shocked at Point C, 
and forms a ring at Point D. (b,c) The variations of the azimuthal and radial 
velocities of the gas along the paths shown in (a). The initial equilibrium 
circular velocity is shown as a dotted line in (b). 

a closed orbit (Point D), forming a nuclear ring at that radius. 
In other words, the centrifugal barrier inhibits the inflowing 
gas from moving further in. Regardless of the BH mass, this 
happens at R ~ (0.8- 1.2) kpc in models with c s = 5 km s -1 
and R ~ 0.4 - 0.6 kpc in models with c s = 20 km s -1 (Fig. 
0]3,c). The facts that the off-axis shocks penetrate the ILR 
in bh7/bh8 models, the outer ILR in bhO models, and that the 
ring positions are almost independent of Mbh when the rings 
are beginning to form suggest that the ring formation is un- 
likely to be governed by resonances. For models with low 
sound speed, the shape of nuclear rings is similar to an x<i- 
orbit. Clearly, the presence of the nuclear ring prevents the 
shocked gas from inf ailing directly to the nucleus. 

Since the off-axis shocks are curved, Crocco's theorem 
ensures that vorticity can be generated at the shock fronts. 
Figure [6] plots snapshots of the potential vorticity £ = | V x 
u + 2^|/E relative to the initial value £o near the shocks at 
t = 90, 100, 110, and 120 Myr of the standard model. At 
t = 90 Myr, £/£ is largest along the shocks. Vorticity pro- 
duced at the shocks is advected with the background flows 
and enters the shock fronts at the opposite side after a half 
revolution. Vorticity grows secularly with time by succes- 
sive passages across the shocks. When vorticity achieves 



substantial amplitudes, it causes the shock fronts to wiggle 
and fragment into small clumps with high vorticity (see Fig. 
Hfc). The process of clump formation along the shocks in our 
models bears remarkable res emblance to the wiggle instabil- 
ity of spiral shocks fo und by IWada & Kodal (|2004|) (see also 
iKim & Os triker 2006). These clumps are carried radially in- 
ward and add to the nuclear ring, making the latter fairly in- 
homogeneous (Fig.|4]l,e). 

The off-axis shocks shown in Figure [5] are not stationary 
largely because the bar potential is not fully turned on yet. As 
the strength of the bar potential keeps increasing, they become 
stronger and move slowly toward the bar major axis. Gas that 
is added to the ring from the off-axis shocks has increasingly 
lower angular momentum with time, causing the ring to shrink 
in radius with time. After the bar potential is fully turned on, 
the off-axis shocks become gradually weaker as the amount 
of gas in the mid-bar regions lost to the ring increases with 
time. At the same time, orbital phase mixing and frequent 
clump collisions in the ring make the latter rounder and align 
its semimajor axis in the direction perpendicular to the bar 
major axis. Note that the rings are always attached to the 
inner end of the off-axis shocks. 

At t = 300 Myr, Model cs05bh0 reaches a quasi-steady state 
in the sense that temporal changes in the overall flow pat- 
tern are very slow, and the locations of the shocks and rings 
do not vary much with time. Figure |4f overplots some of 
the x\ -orbits (dotted curves aligned vertically) and X2-orbits 
(solid curves aligned horizontally), showing that the shape of 
the nuclear ring matches well with an X2 -orbit, while the off- 
axis shocks closely follow one of the x\ -orbits over the whole 
length of the shocks. This is because when c s = 5 km s -1 the 
impact of thermal pressure on the gas orbits is much smaller 
than that of the gravitational and centrifugal forces, so that 
pure orbit theory (neglecting pressure forces) is a good de- 
scription. When c s £ 15 km s -1 , however, thermal pressure 
gradients strongly affect gas orbits, and thus the morphology 
of substructures in the central regions are modified, as we will 
discuss below. 



3.2. Off-axis Shocks 

Even if the gravitational potential is the same, the flow mor- 
phology and velocity fields differ considerably depending on 
the sound speed. Figure shows instantaneous streamlines 
plotted over the logarithm of the density distribution in Mod- 
els cs05bh0 and cs20bh0 at t = 300 Myr. Red curves de- 
note the streamlines that go through the off-axis shocks, while 
those enveloping the off-axis shocks are represented by green 
curves. In all models, the off-axis shocks are almost parallel 
to x\ -orbits. They start from the bar major axis at the outer 
ends, offset toward downstream in the mid-bar regions, and 
connect to the nuclear rings at the inner ends. The mean off- 
set of off-axis shocks from the bar major axis is larger for 
models with smaller c s . 

The outer end regions of the off-axis shocks have com- 
plicated density structures including t he "4/1 -spiral shocks" 
mark ed with a red arrow in Figure (Englmaier & Gerhard 
1997) and the enhanced density ridges (a w hite arrow) termed 
"smudges" bylPatsis & A thanassoula ( 2000). As discussed by 
Englmaier & Gerhard (1997j), the 4/1 -spiral shocks are pro- 
duced by collisions of gas m oving on x\ -orbits with that o n 
the 4/1 -resonant family (e.g., Contopoulo s & Grosb ol 1980). 
The gas loses angular momentum at the shocks and subse- 
quently switches to lower orbits. As the streamlines in green 
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Fig. 6. — Snapshots of the potential vorticity £=|Vxu + 2Q|/E normalized by the initial value £o in Model cs05bh0. Only the regions with -2 kpc < x < 
-0.5 kpc and -1 kpc < y < 2 kpc around the off-axis shocks are shown. Colorba r show s log(£/£o)- This vortex-generating instability of curved shocks appears 
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Fig. 7. — Logarithm of the density distribution overlaid with instantaneous streamlines in Models cs05bh0 and cs20bh0 at t = 300 Myr. The red lines represent 
streamlines that meet the off-axis shocks, while the green lines are for those that go around the shocks. The red and white arrows in (a) mark the 4/1 -spiral shocks 
and "smudge", respectively. The short white line segment in each panel indicates a slit along which density and velocity are measured in Figure[8] 



display, in models with c s = 5 km s -1 , the 4/1 -spiral shocks 
are quite strong and spatially extended, so that the orbits after 
the shocks become relatively radial and converge at the op- 
posite side of the bar, building a smudge after about a half 
revolution. Collisions of streams off the 4/1 -spiral shock and 
the smudge on the same side of the bar funnel the gas at the 
intersections to an x\ -orbit, which are the starting points of 
the off-axis shocks. When c s = 20 km s -1 , on the other hand, 
the 4/1 -spiral shocks are short and weak, and the streamlines 
off the shocks diverge, so that a dense ridge does not form. 
Since the gas becomes less compressible with increasing c s , 



steady off-axis shocks in models with large c s can be sup- 
ported only in inner regions where the bar perturbations are 
sufficiently strong. This explains why the mean offset of the 
off-axis shocks from the bar major axis beco mes smaller as c s 
increases (e.g., Englmaier & Gerhard 1997). With weak 4/1- 
spiral shocks and no smudge, the gas in the bar-end regions in 
model with c s = 20 km s -1 is comparatively unsteady, some- 
times generating small dense blobs that move inward along 
the off-axis shocks. 

To quantify the shock properties, we place a slit in each 
model, indicated by the short white lines in Figure [TJ The 



Nuclear Structures in Barred Galaxies 



TABLE 2 
Properties of Off- axis Shocks 



Model 




V / V 
^max/ ^0 


A A 


A A 
JVi || 


hh 


/ 

dv\\ /as 


Oimax 




(kpc) 








(cleg) 


(10 3 kms -1 kpc -1 ) 




cs05bh0 


0.98 


5.7 


12.1 


10.9 


47.7 


2.7 


7.2 


cslObhO 


0.72 


3.6 


6.0 


4.9 


50.8 


1.8 


3.5 


csl5bh0 


0.55 


5.7 


5.7 


1.0 


80.4 


2.3 


2.6 


cs20bh0 


0.45 


5.0 


3.5 


1.2 


71.3 


1.8 


1.6 


cs05bh7 


0.93 


5.8 


15.4 


6.3 


67.9 


1.8 


7.4 


csl0bh7 


0.63 


9.8 


8.1 


-1.3 


-81.4 


1.6 


4.7 


csl5bh7 


0.47 


9.0 


6.6 


-1.2 


-79.8 


1.7 


3.0 


cs20bh7 


0.31 


5.7 


2.3 


1.3 


59.4 


1.0 


1.9 


cs05bh8 


0.96 


6.1 


16.7 


7.2 


66.7 


2.0 


7.5 


csl0bh8 


0.74 


7.4 


9.6 


-0.1 


-89.5 


2.3 


6.3 


csl5bh8 


0.43 


6.8 


7.2 


-0.8 


-83.3 


1.4 


3.0 


cs20bh8 


0.32 


21.6 


5.9 


-2.6 


-66.6 


1.4 


2.9 



Note. — is the shock position along the slit; S max is the maximum density 
attained immediately after s^; A4± and .Mm are the Mach numbers of the incident 
flow perpendicular and parallel to the shock, respectively; 4h is the inclination angle 
of the incident flow with relative to the shock; dv\\/ds is the velocity shear in the 
postshock region; a max is the maximum value (occurring at the shock front) of the 
compression factor a = -(V • \)AR/c s . 




2.0 



Fig. 8. — Profiles of surface density E, velocity perpendicular and vn 
parallel to the off-axis shock, and compression factor a = -(V • \)AR/c s 
along the slit in bhO models with differing c s at t = 300 Myr. The position of 
the slit is shown in Figure Q 



slit starts from (x,y) = (0,-1.5 kpc) and runs perpendicular 
to the local segment of the off-axis shocks, roughly at R ~ 
(1.5-1.8) kpc. Figure [8] plots the profiles along the slit of 
surface density, velocities, and the compression factor 



a = -(V -\)AR/c s . 



(3) 



the last of which can be used as an effective measure of 
the s hock strength (e.g., Maci ejewskil l2004bt iThakur etaLl 
20090. The gas is flowing from left to right in the increas- 
ing direction of s, where s measures the distance along the 
slit from the starting point. Table [2] gives the shock proper- 
ties for all models at t = 300 Myr: s S h is the position of the 
off-axis shocks along the slit, S max is the peak density after 
^sh, M-± =v±/c s and M\\ = v\\/c s are the Mach numbers of 
the incident flows in the directions perpendicular and parallel 
to the shocks, respectively, 4h = tan _1 (A^±/A^||) is the in- 
clination angle of the preshock velocity relative to the shock 
front, dv\\/ds quantifies the velocity shear in the postshock re- 
gion, and a max is the maximum value of the compression fac- 
tor occurring at the shock front. It is apparent that the off-axis 
shocks tend to move toward the bar major axis with increasing 
c s , while there is no clear dependence of on the BH mass. 
The compression factor at the shock is insensitive to M B h and 
scales roughly with c s as a max ~ 1 .l(c s /5 km s -1 ) 92 . 

Naively, one would expect that the off-axis shocks become 
weaker as the sound speed increases, since the density jump 
in planer isothermal shocks is proportional to -Mj. How- 
ever, this is not the case, as Figure [8] and Table [2] demon- 
strate. There are several reasons for this. Firstly, it is the 
velocity component normal to the shock front v±_ that deter- 
mines the shock jump conditions, and because the inclination 
angle of the streamlines which enter the shock varies with lo- 
cation and with c s , v±_ varies in a complicated fashion. For 
example, Figure [8] shows that for off-axis shocks formed at 
R ~ (1.5-1.8) kpc, Model csl5bh0 with c s = 15 km s" 1 has 
the largest peak density as well as the largest = 85 km s -1 
and z' S h = 80°. On the other hand, Model cslObhO has the 
smallest v± = 60 km s -1 (with / S h = 51°) and thus the lowest 
density enhancement. Since the sound speed is lower, Model 
cs05bh0 with = 60 km s -1 produces S max comparable to 
that in Model csl5bh0. Secondly, we note that the Rankin- 
Hugoniot jump conditions for stationary planar shocks are not 
applicable to the curved and two-dimensional off-axis shocks 
formed in our simulations. As Figure [7] displays, the flows 
are fully two-dimensional in the sense that streamlines di- 

4 For planar isothermal shocks in steady state, a = A4± - M.~^ at the shock 
discontinuities. 
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Fig. 9. — Effects of sound speed and BH mass on the distribution of gas surface density, shown in logarithmic scale, in the central regions of all models at 
t = 300 Myr. The nuclear rings are narrow when c s = 5 km s _1 , while they spread out as c s increases. 



verge before the shocks, and converge after the shock with 
the radial inflow coming from the regions near the end of 
the bar. The fact that the compression factor a measured at 
the shock front is smaller than M±-M~2 expected from 
planer isothermal shocks also indicates that the shocks are 
two dimensional. Finally, the density and velocity fluctua- 
tions generated by the vortex-generating instability are im- 
portant around the off-axis shocks especially for models with 
small c s , so that the flows are not strictly stationary. Note 

5 Negative values of v_l right after the shocks in model cslObhO shown in 



that the shocked gas has strong velocity shear, amounting to 
dv\\/ds ~ (1 -3) x 10 3 km s _1 kpc -1 , which is about 10 times 
larger than the velocity shear arising from galaxy rotation in 
the solar neighborhood. Such strong shear can stabilize the 
high-density, off-axis shocks against self-gravity. 

3.3. Nuclear Rings 

Gas that loses angular momentum at the off-axis shocks 
flows radially inwards and forms a nuclear ring in the cen- 
tral regions. Figure [9] shows diverse morphological features 

Fig-EJare due to vortices produced by the instability. 
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Fig. 10. — Instantaneous streamlines (red solid lines) overlaid on the logarithm of the density distribution for models with (a) c s =5 km s and no BH, (b) 



c s = 20 km s -1 and no BH, (c) c s = 20 km s~ 
panels represent X2 -orbits. 



and M B h = 4 x 10 7 M , and (d) c s = 20 km s -1 and M B n = 4 x 10 8 M , at t = 300 Myr. The dotted curves in all 



TABLE 3 
Properties of Nuclear Rings 



Model 


R m (kpc) 


Rout (kpc) 


7? ring (kpc) 


(E) ma x/Eo 


^ring/So 


cs05bh0 


0.63 


1.22 


0.92 


17.2 


11.9 


cs05bh0t 


0.62 


1.17 


0.89 


18.4 


13.1 


cslObhO 


0.45 


0.90 


0.68 


30.2 


17.9 


csl5bh0 


0.30 


0.77 


0.54 


37.5 


17.6 


cs20bh0 


0.23 


0.71 


0.49 


36.8 


16.9 


cs20bh0t 


0.20 


0.73 


0.48 


34.7 


16.1 


cs05bh7 


0.66 


1.17 


0.90 


18.9 


13.9 


csl0bh7 


0.46 


0.88 


0.67 


30.3 


18.9 


csl5bh7 


0.28 


0.80 


0.52 


33.9 


16.9 


cs20bh7 


0.23 


0.67 


0.47 


41.0 


19.5 


cs20bh7t 


0.22 


0.66 


0.46 


43.4 


19.9 


cs05bh8 


0.69 


1.25 


0.94 


18.1 


12.7 


csl0bh8 


0.47 


0.89 


0.67 


30.9 


19.4 



csl5bh8 
cs20bh8 



Note. — R m and R ou t are the inner and outer radii of the ring defined 
by the positions where (E) = (E) max /5, with (S) max being the maximum 
density; R r - mg is the mass-weighted ring radius; E r i ng is the mean density of 
the ring. 



produced in the regions with \x\,\y\ < 1 kpc for all models at 
t = 300 Myr. Figure [TO] overplots instantaneous streamlines 
for a few selected models. Some models (with low c s ) have a 



nuclear ring together with inner spiral structures, some mod- 
els (with high c s and small M B h) have a ring with no spirals, 
while others (with high c s and large M B h) possess only nu- 
clear spirals without an appreciable ring. 

When c s = 5 km s -1 , the nuclear rings are quite narrow and 
clearly decoupled from the nuclear spirals. Even though the 
ring has a large density, the low sound speed makes the ef- 
fect of thermal pressure on the gas orbits insignificant. The 
gas around the ring in Model cs05bh0 thus follows X2-orbits 
fairly well and the shape of the ring does not deviate con- 
siderably from X2-orbits (Fig. fTOh). When c s = 20 km s -1 , 
on the other hand, the pressure force in the central regions 
becomes important and affects the shape of the gas stream- 
lines. Even the inflowing gas that arrives at the contact points 
between the off-axis shocks and the nuclear ring takes very 
different orbits depending on its location. Some gas at the 
outer parts of the contact points is pushed out by the pres- 
sure gradient and follows trajectories that are much rounder 
than X2 -orbits, while the gas in the inner parts is forced to take 
inner highly eccentric orbits (Fig. [T0b). Consequently, the 
gas in the central regions in Models cs20bh0 spreads out spa- 
tially and forms a ring that is more circular and broader than 
in Model cs05bh0. Since the presence of a central BH in- 
creases the initial angular momentum of the gas in the central 
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Fig. 11. — Radial distribution of the gas surface density (E) averaged 
both azimuthally and temporally over t = 300 - 500 Myr for models with (a) 
M BH (0) = 0, (b) Mbh(0) = 4 x 10 7 M , and (c) M BH (0) = 8 x 10 7 M . The 
dotted lines are the results of the models in which Mbh varies with time. 
The locations of Rm&x and R m [ n where the Q-k/2 curve attains local maxi- 
mum and minimum and the relevant ILRs are indicated as arrows along the 
abscissa. In (c), the dashed lines correspond to the cases with c s = 15 or 
20 km s _1 for which the density at R < 0.1 kpc is dominated by nuclear spi- 
rals rather than rings. 

regions, the pressure effect becomes less important as M B h 
increases. Figure [TO1 shows that the pressure distortion of x^- 
orbits is still significant for M B h = 4 x 10 7 M , while the gas 
orbits in the very central parts at R £ 0.2 kpc remain almost 
intact when M B h = 4 x 10 8 M . In Model cs20bh8, some 
gas at7?~0.5-0.8 kpc temporarily moves radially outward 
due to the radial pressure gradient built up by the background 
gas and is subsequently swept inward by other gas flowing in 
along the off-axis shocks. 

Figure [TT] plots the radial distribution of the averaged gas 
surface density (E), averaged both azimuthally and tempo- 
rally over t = 300-500 Myr for all models. The locations 
of the ILRs as well as R mSLX and R m [ n corresponding to the 
local maximum and minimum of the Q - kJ2 curve are in- 
dicated by arrows on the abscissa. Table [3] gives the inner 
and outer radii, R[ n and R 0VLt , of the ring, the mass-weighted 

ring radius R ring = J^ ut (E)RdR/ J^ ut (Y,)dR, the peak density 

(E) max , and the mean density E ring = J^ out (Z)dR/(R 0Ut -R in ) 
of the ring in each model. Here, R{ n and R 0VLt are defined as the 
radii where (E) = (E) max /5. Note that Models csl5bh8 and 



cs20bh8 do not harbor a well-defined nuclear ring, as will be 
discussed in the next section. 

All rings that form are located within 7?oilr if there are 
two ILRs or 7? ILR if there is a single ILR, indicating that the 
formation of a nuclear ring does not require the presence of 
two ILRs. However, there is in general no direct connec- 
tion between the ring positions and 7?oilr or 7? ILR . When 
c s = 5 km s -1 , the rings are all located at R ~ 0.6- 1.2 kpc, 
independent of Mbh- The mass- weighted radius is Rn ng ~ 
(0.90-0.92) kpc, indicating that the ring position is not gov- 
erned by the shape of the ft - k/2 curve. When c s = 10 km s -1 , 
the ring radius decreases to R Y { ng ~ (0.67-0.68) kpc, again 
insensitive to the BH mass, consistent with the tendency for 
the off-axis shocks to move closer to the bar major axis as 
c s increases. Rings with c s = 10 km s -1 have a larger sur- 
face density than those with c s = 5 km s -1 , corresponding ap- 
proximately to a constant ring mass (i.e., E r i ng ex R~^ g ). As 
c s increases further, high thermal pressure provides strong 
perturbations for X2 -orbits and tend to spread out the gas in 
the central parts, resulting in a broad distribution of (E) at 
R <, 0.5 kpc. When the BH is not massive enough, these 
perturbations wipe out coherent, weak spiral structures that 
formed earlier in the nuclear regions. 

Because the presence of a central BH dominates the poten- 
tial only in the central region, the allowance for the growth 
of the BH due to mass accretion does not make significant 
difference in the regions outside the ring. Even inside the 
ring, Figure [TT] and Table [3] show that the changes in the ring 
size 7? r i n g and the ring density E r i ng caused by the temporal 
change of Mbh are less than 4% and 10%, respectively. We 
will show below that BH growth in our simulations does not 
significantly affect the properties of nuclear spirals and mass 
inflow rates, as well. 

4. NUCLEAR SPIRALS 

High-resolution observations of barred galaxies reveal that 
some contain nucle ar spirals inside a nuclear ring (e.g.. 
Marti ni et al.ll2003allbl: iPrieto et alJl2005 Ivan de Yen & Fathll 
2010). As Figures [9] displays, some of our models also have 
spiral structures in their nuclear regions that persist for long 
periods of time. Other models also have nuclear spirals at 
early time but they are destroyed as a result of interactions 
with nuclear rings. In this section, we describe the formation 
and shape of nuclear spirals in detail. Figure [12] schematically 
summarizes how the type of nuclear spirals changes with time 
and how long they survive in each model. Table [4] lists the 
properties of nuclear spirals measured at t = 500 Myr for the 
models that possess long-lasting spirals: R{ n and R out denote 
the radii of the inner and outer ends of the nuclear spirals, re- 
spectively; E avg and E pea k are the azimuthally-averaged and 
peak surface densities, respectively, at R = R a = 0.05 kpc for 
Models csl5bh8 and cs20bh8 and R a = y/RmR ut for the other 
models; i p is the pitch angle of the spirals at R = R a . A nega- 
tive (positive) value of i p indicates leading (trailing) spirals. 

4.1. Models Without A Black Hole 

To study the spiral features, it is convenient to show the 
logarithm of the gas surface density in logarithmic polar co- 
ordinates. Figure [13] plots snapshots of gas surface density of 
Models cs05bh0 and cslObhO on the (j)-logR plane. Any co- 
herent features with a positive (negative) slope on this plane 
are leading (trailing) waves. Only the regions with R < 1 kpc 
are shown. Two horizontal lines mark Rulr and R max where 
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Fig. 12. — Schematic illustration of the types of nuclear spirals found in our simulations, and their duration. 
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Fig. 13. — Snapshots of the logarithm of the gas surface density on the $-\ogR plane for Models cs05bh0 (top row) and cslObhO (bottom row). These models 
do no have a central BH. Only the regions with R < 1 kpc are shown. The locations of Rulr and R m ax are indicated by two horizontal lines. 
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TABLE 4 
Properties of Nuclear Spirals 



Model 


R in (kpc) 


Rout (kpc) 


^avg/^0 


^peak/^avg 


i P (deg) 


cs05bh0 


0.13 


0.50 


0.41 


7.92 


-33.8 


cs05bh0t 


0.20 


0.45 


0.56 


9.84 


-35.4 


cs05bh7 


0.02 


0.45 


0.92 


1.95 


8.5 


csl0bh7 


0.05 


0.42 


0.16 


3.86 


55.1 


cs05bh8 


0.02 


0.55 


1.00 


2.07 


3.5 


csl0bh8 


0.02 


0.40 


2.28 


1.91 


5.8 


csl5bh8 


0.02 




8.09 


1.67 


6.2 


cs20bh8 


0.02 




170.0 


1.74 


7.8 



Note. — R in and R out are the inner and outer ends; Savg and Sp ea k 
are the mean and peak densities at R a = 0.05 kpc for Models csl5bh8 and 
cs20bh8 and R a = CKin^out) 1 / 2 for the other models; i p is the pitch angle 
atR = R a . 



the Q — k/2 curve is a locally maximum. At early time, the 
non-axisymmetric bar potential induces weak m = 2 pertur- 
bations in the central regions. Perturbed gas elements in 
the galactic plane follow slightly elliptical orbits, which are 
closed in a frame rotating at ft - n/2 and thus precess at a rate 
Q — k/2 near the ILRs whe n seen in a stationar y bar frame. 
As succinctly depicted in Buta & Combes) (119961) . due to col- 
lisional dissipation of gas kinetic energy occurring on con- 
verging orbits, the gas forms spiral structures whose shape de- 
pends critically on the sign of d(Vt -n/2)/dR such that spirals 
are leading (trailing) where d(Q - n/2)/dR is positive (nega- 
tive). Figure [13] indeed shows that when t = 50 Myr the per- 
turbed density is leading at R < R max and trailing at R > 7? max , 
although the trailing features are soon overwhelmed by the 
nuclear ring. Located away from the nuclear ring, however, 
the inner leading waves are able to grow with time and eventu- 
ally develop into shock waves at t ~ 200 Myr. These nuclear 
spirals are short, extending over R ~ 0.13-0.50 kpc, quite 
open with a pitch angle of i p = -30°, and almost completely 
detached from the nuclear ring. 

Figure [14] plots the azimuthal distributions of surface den- 
sity and velocities of the nuclear spirals at R = 0.25 kpc 
in Models cs05bh0 (solid lines) and cs05bh0t (dotted lines) 
when t = 500 Myr. Over the course of the orbits, the 
changes of the radial and azimuthal velocities associated with 
the spirals amount to ~ 100 km s -1 , which is indeed large 
enough to induce shocks. The peak densities occurring at 
~ 100° , 280° for Model cs05bh0 correspond to shock fronts 
with a compression factor of a ~ 3.4. The density bumps at 
(j) ~ 135°, 315° are produced by waves launched from the in- 
ner boundary. In Model cs05bh0t, the BH mass is increased 
to M B h ~ 10 5 M due to mass inflow, which supports slightly 
stronger, more leading spiral shocks than in Model cs05bh0 
with no BH. Despite continual perturbations by traveling trail- 
ing waves propagating from the inner boundary, the nuclear 
spirals in these models last until the end of the run. That 
leading nuclear spirals are persistent when the sound speed 
is small i s consistent with the r esults of SPH simulations re- 
ported bv lAnn & Thakud (12005b . 

The usual WKB dispersion relation for tightly-wound, 
linear-amplitude waves in a non-self-gravitating medium 
reads 

(uj-mn) 2 = k 2 c 2 s +K 2 (4) 

where uj is the wave frequency and k and m are the radial an d 
azimuthal wavenumbers (e.g., lGoldreich & Trema ine1ll979|) . 
For m = 2 waves corotating with a bar (i.e., uj = 2Qb), equation 




-loo t , , , i , , , i 

180 360 

( d eg) 

FIG. 14. — Azimuthal profiles of surface density (top) and velocities (bot- 
tom) of the nuclear spirals at R = 0.25 kpc in Models cs05bh0 (solid) and 
cs05bh0t (dotted) when t = 500 Myr. The spirals at this radius are shocks 
with the compression factor of a ~ 3.4. 

© becomes 



k = ± — + k/2 - n b )(n - n/2 - Q b ), (5) 

c s 

(e.g.. lEnglmaier & Shlosmanl 120001: lMacieiewskill2004ah . in- 
dicating that nuclear spirals corotating with the bar can ex- 
ist only in the regions between R\jlr and 7?oilr when there 
are two ILRfl However, the top row in Figure [13] reveals 
that the nuclear spirals in Model cs05bh0 extend slightly in- 
ward of Rhlr- This seemingly contradicting result is due 
to nonlinear effects which are not captured by WKB the- 
ory. The velocity perturbations associated with these spi- 
rals are so large that fluid elements just outside Rulr can 
move across the inner ILR over the course of their epicy- 
cle orbits, providing perturbations for the gas at R < Rulr 
that responds passively. In Model cs05bh0, the radial veloc- 
ity perturbation is Avr = 52 km s" 1 . Since the epicycle fre- 
quency is k = 1 100 km s -1 kpc -1 at R = Rulr, the correspond- 
ing radial amplitude of the epicycle orbits is estimated to be 
AR = Avr/k = 0.05 kpc, which is in good agreement with the 
radial extent of the nuclear spirals inward of Rulr. 

Models without a BH and with c s > 10 km s -1 do form nu- 
clear spirals at early time, but they are all transient, lasting less 
than 200 Myr. The bottom row of Figure[9] shows how nuclear 
spirals are destroyed in Model cslObhO. The nuclear ring in 
this model is not only located inside R max but also has large 
thermal pressure, continuously generating sonic perturbations 
that propagate radially inward. Because of the background 
shear, the perturbations are preferentially in the form of trail- 
ing waves which interact destructively with the leading spirals 
that formed at R < R mSLX , destroying the latter. The destruction 
of nuclear spirals happens at t = 185 Myr for Model cslObhO. 

6 Fig.[l3]shows that there are low-amplitude waves propagating relative to 
the bar inside Rulr. Such waves can exist inside the inner ILR as long as 
they satisfy equation in the WKB limit. 
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This occurs earlier as c s increases, since disturbing pressure 
perturbations are correspondingly stronger. 

4.2. Models With M BH = 4 x 10 7 M 

Since the presence of a central BH greatly changes the 
Q — k/2 curve in the central regions, it is interesting to ex- 
plore how the morphologies of nuclear spirals depend on the 
BH mass. In bh7 and bh8 models with a single ILR, sta- 
tionary waves in the bar frame can exist inside T^lr all the 
way down to the center. Figure [15] plots snapshots of gas sur- 
face density for Models cs05bh7 and csl0bh7 on the loga- 
rithmic polar plane. The positions of R max and R m [ n are in- 
dicated by the horizontal lines. As expected, the overdense 
perturbations produced by orbit crowding at early time it = 50 
Myr) have leading configurations at R m m <R< Rmax and trail- 
ing configurations at R < R m ] n or R > R max . The overdense 
regions at R > R max are subsequently wiped out as the nu- 
clear ring forms, while those at R < R max grow into trail- 
ing nuclear spirals. In Model cs05bh7, the leading spirals 
at R m m < R < Rmax also grow slightly until t ~ 150 Myr to 
temporarily form "inner- trailing and outer-leading" structures 
represented by the double cross-hatching in Figure [12] As 
trailing perturbations from both the inner trailing spirals and 
the outer nuclear ring propagate and interfere with the leading 
spirals, it becomes increasingly difficult to identify coherent 
spiral structures at < R < R m2LX - On the other hand, the 
trailing spirals at R < R m [ n keep growing until t ~ 230 Myr 
after which their amplitude of £ pe ak/£avg ~ 2.0 remains more 
or less constant. They are approximately logarithmic in shape 
with a pitch angle of i p = 8.5°. Unlike in Model cs05bh0 
where leading spirals are actually shocks, the trailing nuclear 
spirals in Models cs05bh7 are relatively weak (with the max- 
imum compression factor of a max ~ 0.28 at t = 500 Myr), and 
never develop into shocks. 

In Model csl0bh7 (bottom row of Fig. [15]), the nuclear spi- 
rals have larger k and thus are more open than those in the 
c s = 5 km s -1 counterparts (see, e.g., eq. 0). Since the nu- 
clear ring in this model forms at Rn ng ^0.5 kpc, it can di- 
rectly destroy the leading spiral at R m [ n <R< R max , and feeds 
the inner trailing spirals by supplying trailing perturbations. 
Thus, the inner spirals grow both in strength and spatial ex- 
tent to make contact with the nuclear ring at t = 210 Myr. At 
this time, all parts of the nuclear spirals turn to shocks with 
the maximum density of £ max /Eo = 3.8 and the correspond- 
ing compression factor of a = 0.8 at R = 0.1 kpc. Gas passing 
through the spiral shocks loses angular momentum, increas- 
ing the mass inflow rate at the inner boundary. As the amount 
of gas lost in the nuclear regions increases, the density of the 
nuclear spirals decreases with time, but they remain as shocks 
with the compression factor of a ~ 1.3 until the end of the 
run. 

In Models csl5bh7 and cs20bh7, nuclear spirals start out 
with inner- trailing and outer-leading shapes, and evolve into 
trailing-only configurations, as in the other bh7 models. How- 
ever, the highly eccentric orbits of the gas affected by ther- 
mal pressure dismantle the inner spiral structures almost com- 
pletely in these models: thermal perturbations are so strong 
that a central BH with Mbh = 4 x 10 7 M cannot enforce cir- 
cular orbits in the very nuclear regions (see Fig.fTOb). 

4.3. Models With M BH = 4 x 10 8 M 

Since the ft - k/2 curve decreases monotonically with R < 
1 kpc in bh8 models, nuclear spirals, if they exist, are all trail- 
ing as evident in Figure [16] For Model cs05bh8 (top row 



of Fig. [16]), the spirals evolve almost independently of, and 
are well separated from, the nuclear ring. At t = 500 Myr, 
they have a very small pitch angle i p = 3.5° corresponding to 
kR = 2 cot i p = 33 owing to the strong background shear. They 
are also very weak, with an amplitude of S pea k/E avg = 2.1 at 
R « 0. 1 kpc. When c s = 10 km s -1 , the pitch angle is increased 
to i p = 5.8° and the nuclear spirals extend outward all the way 
to the nuclear ring. Except near the contact points, the spirals 
are still decoupled from the ring. With quite a large value of 
kR = 20, the component of the rotational velocity perpendicu- 
lar to the spirals is not large enough to induce shocks. 

We have seen earlier that the large thermal pressure in the 
rings of models with c s > 15 km s -1 provides strong pertur- 
bations that make the gas orbits near the galaxy center highly 
eccentric, destroying nuclear spirals in bhO and bh7 models. 
In bh8 models, however, the situation is quite different since 
a central BH with M B h = 4 x 10 8 M dominates the potential, 
keeping the orbits almost circular in the very central regions. 
Even with a large thermal pressure, the gas orbits there can- 
not be very eccentric, so that the nuclear spirals are protected 
from disruptive pressure perturbations. On the other hand, the 
ring material is quite distributed because of the large pressure 
gradients, feeding a trailing nuclear spiral that grows strongly 
in Models cs 1 5bh8 and cs20bh8 . 

At t = 120 Myr, the spirals in Models csl5bh8 and cs20bh8 
turn into shocks and touch the densest parts of the ring located 
at R ~ 0.4 kpc. Because of the larger pitch angles, the shocks 
are stronger in the outer parts; the portions at R £ 0. 1 kpc are 
weak shocks with density jumps of only ~ 2. Similarly, the 
shocks in Model cs20bh8 are stronger than in Model csl5bh8 
since the former has more open spirals and a denser ring. In 
both models, the shocks near the ring are so strong that even 
the gas constituting the ring suffers from a significant loss of 
angular momentum at the intersecting points. At t = 150 Myr, 
the ring material is essentially dissected by the trailing spiral 
shocks and gradually moves toward the center. As the ring 
material continues to flow in, the spirals appear as a direct 
continuation of the off-axis shocks (t = 350 Myr), consistent 
with the results of Maciejewski ( 20042). 

Figure [TTlplots the radial distributions of the maximum and 
mean densities as well as the maximum compression factor 
in the inner 1 kpc regions of Models csl5bh8 and cs20bh8 at 
t = 500 Myr. The inflow of the ring material in Model cs20bh8 
is quite strong that the gas is collected in the nuclear regions 
with R ~ 0.02-0.1 kpc; the amount of the gas that goes out 
through the inner boundary is much smaller than that comes 
in. With weaker shocks and thus less angular momentum loss, 
on the other hand, the destroyed ring gas in Model csl5bh8 is 
still mostly at R > 0.1 kpc at this time. Note that the den- 
sity enhancement at R = 0.05 kpc resulting from the dissolu- 
tion of the ring is ~ 8 So and ~ 170£o m Models csl5bh8 
and cs20bh8, respectively: the mean contrast of the spirals 
that are weak shocks with the compression factor a ~ 0.3 is 
Speak/S a vg ~ 1 .6 - 1 .8 at R = 0.05 kpc in both models. This 
increase of the gas surface density in the central regions is the 
primary reason for enhanced mass inflow rates at late time in 
bh8 models with large c s . 

5. MASS INFLOW RATES 

Galactic bars are considered to be a promising means of 
transporting gas to the centers of galaxies to fuel supermas- 
sive BHs and produce AGNs. Since our numerical models 
use a cylindrical grid with a circular boundary, they are ide- 
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Fig. 15. — Snapshots of the logarithm of the gas surface density on the (f)-\ogR plane for Models cs05bh7 (top row) and csl0bh7 (bottom row). These models 
have a central BH with a mass of4x 10 7 M . Only the regions with R < 1 kpc are shown. The locations of R m ax and R m i n are indicated by two horizontal lines. 

ally suited to study how the mass accretion rate depends on 
the gas sound speed and the BH mass. We assume that all 
the gas that crosses the inner boundary located at R = 20 pc 
in our models is accreted to the central BH. In reality, a large 
amount of mass in gas near the BH may change the gas orbits 
by providing pressure and gravitational forces. Therefore, M 
through the inner boundary that we measure is likely to be 
an upper limit to the real mass accretion rate to the BH. In 
addition, M is likely to depend on the inner boundary size es- 
pecially when gas orbits are highly eccentric near the inner 
boundary. Note that in non-self-gravitating, isothermal, and 
unmagnetized systems, M resulting from simulations is lin- 
early proportional to the adopted initial surface density £o; 
all the models presented here take So = 10 M pc~ 2 . 

In general, M would be large if the gas orbits near the cen- 
ter are highly eccentric or radial, while circular orbits would 
make M quite small. This expectation is consistent with Fig- 
ure [18] which plots the temporal evolution of the mass inflow 
rates for all models. The corresponding accreted gas mass 
M acc = / M(t)dt over 500 Myr is given in Table \5\ Clearly, 

M is larger for models with larger c s and no BH, compared small values of M ( < 1O" 3 M yr" 1 ). In bh8 models, the pres- 
to models with a massive BH of M BH = 4 x 10 8 M . When ence of a central BH makes M smaller b Y about an order of 
c s = 5 km s"\ the nuclear spirals are well separated from the magnitude than in bhO models by providing strong axisym- 
nuclear rings, and the departure of the gas orbits from a circu- metric gravitational potential near the center In bh7 models, 

lar shape near the inner boundary is small, resulting in quite the P otential 18 not stron § enou S h t0 ™ larize the eccen " 

trie orbits, giving rise to M only slightly smaller than that in 



Model 


M aC c(M ) 


cs05bh0 


8.9 x 10 4 


cs05bh0t 


1.1 x 10 5 


cslObhO 


8.0 x 10 5 


csl5bh0 


2.7 x 10 6 


cs20bh0 


1.8 x 10 7 


cs20bh0t 


2.3 x 10 7 


cs05bh7 


4.2 x 10 4 


csl0bh7 


1.6 x 10 6 


csl5bh7 


7.4 x 10 6 


cs20bh7 


3.2 x 10 7 


cs20bh7t 


2.9 x 10 7 


cs05bh8 


9.4 x 10 3 


csl0bh8 


4.5 x 10 4 


csl5bh8 


2.0 x 10 5 
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5.1 x 10 6 
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FIG. 16. — Snapshots of the logarithm of the gas surface density on the </>-log 
have a central BH with a mass of 4 x 10 8 Mq . Only the regions with R < 1 kpc 

bhO models. 

Increasing c s enhances M because pressure perturbations 
become stronger and the nuclear ring tends to be located 
closer to the center, both of which strongly affect the gas or- 
bits in the central region. When c s = 10 km s -1 , M is increased 
by about an order of magnitude compared to the cases with 
c s = 5 km s -1 , but is still less than 1CT 2 M yr -1 except for 
a brief time interval around t = 300 Myr in Model csl00bh7 
when the nuclear spirals shock the central gas and cause en- 
hanced inflows. When c s > 15 km s -1 , the gas orbits in bhO 
and bh7 models are quite eccentric near the center, so that 
some gas with highly radial orbits plunges directly into the 
inner boundary, increasing M dramatically compared to the 
lower c s counterparts. The associated gas mass accreted to 
the galaxy center is of order of ~ 10 7 M over 500 Myr, sug- 
gesting that a strong galactic bar like studied in this work can 
be an appealing means for the growth of supermassive BHs 
provided the gaseous medium has a large (e ffective) sound 
speed (cf. Shlosman et al] IT989t see also, e.g., Volonteri 2010 
for review). In bh8 models, on the other hand, the gas orbits in 
the central parts are not perturbed much (since the initial an- 
gular momentum is overwhelmingly large). As the nuclear 
spirals grow into shocks, however, the nuclear gas as well 
as the ring material in these models drift inward slowly, in- 
creasing M over time. In models with Mbh = 4 x 10 8 M , the 
late-time values of M in models with c, = 15 or 20 km s -1 are 



R plane for Models cs05bh8 (top row) and cs20bh8 (bottom row). These models 
are shown. 

larger than 0.0 1 M^ yr" 1 , sufficient to power AGNs in S eyfert 
galaxies (e.g.. lFriedli & Benzlll993tlFabian et al.ll2008l) . 

Finally, we present the mass inflow rates resulting from the 
models in which M B h is varied self-consistently with M. The 
left panels of Figure [19] compare M from Models cs05bh0t 
and cs20bh0t with those from the fixed-MeH counterparts, 
while the bottom panels plot the temporal evolution of the 
BH mass in the former models. In Model cs05bh0t, the total 
gas mass accreted over 500 Myr is ~ 10 5 M , with the corre- 
sponding increment of the equilibrium rotational velocity of 
~ 4 km s -1 at the inner boundary. Since this is three times 
smaller than the initial circular velocity there, the BH mass 
does not affect M much, as Figure [T9l illustrates. In the case 
of Model cs20bh0t, M BH attains - 2 x 10 7 M at t = 500 Myr, 
which is large enough to be dynamically important. When 
c s = 20 km s -1 , however, the gas orbits are highly eccentric and 
M is insensitive to the BH mass as long as Mbh ~ 4 x 10 7 M 
(see Fig. fT8ti). As the rignt panels of Figure [19] show, the in- 
crease of Mbh over 500 Myr in Model cs20bh7t is less than a 
factor of two, corresponding to the equilibrium circular veloc- 
ity 1.3 times the initial value. With an enhanced centrifugal 
barrier, the resulting M becomes gradually smaller than the 
case with fixed M B h, but by less than, on average, ~ 10% in 
t = 300-500 Myr. Therefore, we conclude that the effect of 
BH growth due to gas accretion on the mass inflow rate as 
well as bar substructures is not significant for the models we 
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Fig. 17. — Radial distributions of the maximum and mean densities (top) 
and the maximum compression factor (bottom) in the inner 1 kpc regions of 
Models csl5bh8 and cs20bh8 at t = 500 Myr. Due to strong spiral shocks, the 
ring material in Model cs20bh8 is already moved to the R ~ 0.02-0.1 kpc 
region, while with weaker shocks it still lies at R > 0.1 kpc in Model csl5bh0. 
In both models, the shocks at R ~ 0.05 kpc are quite weak with the compres- 
sion factor of a ~ 0.3. 
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Fig. 1 8 . — Temporal evolution of the mass inflow rates M through the inner 
boundary at R = 20 pc for models with (a) c s =5 km s _1 , (b) c s = 10 km s _1 , 
(c) c s = 15 km s _1 , and (d) c s = 20 km s _1 . The dashed, solid, and dotted 
lines correspond to the cases with M BH = 0, 4 x 10 7 M , and 4 x 10 8 M , 
respectively. 
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Fig. 19. — Temporal evolution of the mass inflow rates (top) and the BH 
mass (bottom) in Models cs05bh0t, cs20bh0t, and cs20bh7t where Mbh is 
allowed to vary with time. The mass inflow rates from the fixed-MBH coun- 
terparts are compared as dotted lines. In all models, the total increase in the 
BH mass over 500 Myr is not large enough to cause significant changes in M. 



have considered. 

6. SUMMARY AND DISCUSSION 
6.1. Summary 

We have presented detailed numerical models that ex- 
plore the formation of substructures produced by the gas 
flow in barred gal axies. Previous models based on particle 
simulations (e.g.,|E nglmaier & Gerhard 19971 : 1 Ann & Th akur 
120051: iThakur et al.1 120091) did not have sufficient resolution 
to resolve nuclear spirals. On the other hand, studies that 
used the grid-based code CMHOG (e.g., iPiner et al.1 1 19951: 
Maciejewski 2004b) unknowingly made mistakes in the force 
evaluation for the bar potential, so that the results needed to 
be recomputed. 

In this paper we have corrected the errors in the origi- 
nal CMHOG code and run high-resolution hydrodynamical 
simulations. To resolve the nuclear regions, we employ a 
logarithmically- spaced cylindrical grid, with a zone size of 
AR < 6 pc at R < 1 kpc where nuclear rings and spirals form. 
We have included the potential from a central BH, and studied 
the flow properties as the mass of the BH Mbh and the sound 
speed c s in the gas are varied. For simplicity, the effects of 
gaseous self-gravity and magnetic fields are not included. The 
main results of the present work are summarized as follows: 

1. Off-axis Shocks - The imposed non-axisymmetric 
bar potential provides gravitational torques to the otherwise 
circular-rotating gas, perturbing its orbit. The perturbed or- 
bits crowd at the downstream sides of the bar major axis and 
produce overdense ridges that eventually develop into off-axis 
shocks. At a quasi-steady state, the off-axis shocks are overall 
almost parallel to x\ -orbits: they start from the bar major axis 
at the outer ends, are gradually displaced downstream as they 
move inward, and connect to the nuclear r^ng at the inner ends. 
While the positions of the off-axis shocks are almost indepen- 
dent of M B h, since the effect of a BH is negligible at large 
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radii, they depend on c s in such a way that the shocks are, on 
average, located closer to the bar major axis as c s increases. 
This is primarily because gas with larger c s should be more 
strongly perturbed to induce shocks, which occurs deeper in 
the potential well and thus re sults in shocks on lower x\ -orbits 
(Engl maier & Gerhardlll997|) . 

The off-axis shocks are in general curved. Flow streamlines 
are complicated near the shocks in that they diverge before 
the shocks and are promptly swept inwards by inflowing gas 
right after the shocks. Therefore, the usual Ranking-Hugoniot 
jump conditions for planar, one-dimensional shocks are not 
applicable to the off-axis shocks. In fact, the shock strength, 
as measured by the peak density E max at ~ (1.5 — 1.8) kpc 
from the center, is £ max /Ho ~ 3 - 6 for models with no BH and 
do not sensitively depend on c s . The compression factor of the 
off-axis shocks are insensitive to the BH mass and depends 
on c s roughly as a max ~ l.l(c s /5 km s -1 ) 92 . The off-axis 
shocks have very strong velocity shear amounting to ~ (1 - 
3) x 10 3 km s -1 kpc -1 . This strong shear may suppress the 
growth of gravitational instability in the high-density off-axis 
shocks when self-gravity is included. 

2. Nuclear Rings - When gas passes through the off-axis 
shocks, it loses angular momentum, flows inward, and forms a 
nuclear ring where the centrifugal force balances the external 
gravity. The ring is attached to the inner ends of the off-axis 
shocks and thus becomes smaller in size as c s increases. The 
rings that form in our models are all located inside the (outer) 
ILR, but this does not imply that the ring formation is related 
to the ILRs. When c s is small, the pressure perturbations on 
the gas orbits are so weak that the rings are quite narrow and 
their shape is well described by x^ -orbits. The mean radius is 
^ring ~ (0.8 - 0.9) kpc when c s = 5 km s -1 and Rn ng ~ (0.5 - 
0.6) kpc when c s = 10 km s" 1 , independent of the BH mass. 
This suggests that the ring position is not determined by the 
Q-k/2 curve, hence by the ILRs. When c s > 15 km s -1 , 
on the other hand, large thermal pressure strongly affects the 
gas orbits in the nuclear rings. For example, some gas near 
the contact points between the ring and the off-axis shocks is 
forced out to follow relatively round orbits, while other gas is 
pulled in radially to make very eccentric orbits. These diverse 
gas orbits near the contact points tend to spread out the ring 
material, making the rings much broader than in models with 
smaller c s . 

3. Nuclear Spirals - Since even the non-axisymmetric 
bar potential is nearly axisymmetric in the central parts, the 
gaseous responses inside a nuclear ring are not as dramatic 
as in off-axis shocks. Nevertheless, non-axisymmetric m = 2 
perturbations are able to grow inside a ring and develop into 
nuclear spirals that persist for a long period of time, provided 
that c s is small or M B h is large. Although all models have 
weak spiral structures at early time, in models with large c s 
they are soon destroyed by eccentric gas orbits as well as per- 
turbations induced by the large pressure in the rings unless 
the BH mass is very large. When M BH = 4 x 10 8 M , the dis- 
ruptive pressure perturbations from the rings cannot penetrate 
the very central parts where the gas has extremely large ini- 
tial angular momentum. In this case, the nuclear spirals are 
protected from the surrounding and thus are long-lived. 

The shape of nuclear spirals is determined by the sign of 
d(Sl -n/2)/dR such that spirals that form in the regions where 
d(Q- n/2)/dR is positive (negative) are l eading (trailing), 
confirming the theoretical expectations of iButa & Co mbes 
(|1996|) . With no BH, only the model with c s = 5 km s -1 has 



persistent leading spirals in the regions where the Q-k/2 
curve is an increasing function of R. The leading spirals 
in this model are quite strong and develop into shocks with 
the peak density X pe ak/£avg ~ 7.9 and the compression fac- 
tor a ~ 3.4 at R = 0.25 kpc at the end of the run. Models 
with Mbh = 4 x 10 7 M initially have hybrid features com- 
prising of trailing spirals at R < R m m and leading spirals at 
^min <R< Rmax, where R m [ n and R max refer to the radii where 
Q-k/2 curve attains a local minimum and maximum, respec- 
tively. When c s < 10 km s -1 , however, the leading parts in 
the hybrid spirals are destroyed by the trailing pressure waves 
launched by the ring, leaving only the weak trailing spirals 
behind. When c s > 15 km s" 1 , both leading and trailing spi- 
rals are destructed completely by the pressure perturbations. 
In models with M BH = 4 x 10 8 M , d(Q-K,/2)/dR < in the 
whole nuclear regions, so that nuclear rings are always trail- 
ing. When c s < 10 km s -1 , the spirals well separated from the 
ring are weak and tightly wound. When c s > 15 km s -1 , on 
the other hand, the trailing spirals are fed by the gas in the 
nuclear ring and grow into shocks, the outer ends of which 
join the inner ends of the off-axis shocks smoothly. Although 
the nuclear spirals in these models have only modest den- 
sity contrasts of S P eak/S av g ~ 1.7 at R = 0.05 kpc, the back- 
ground density resulting from the inflows of the ring material 
is greatly enhanced (by more than two orders of magnitude 
when c s = 20 km s -1 ). 

4. Mass Inflow Rates - Although gas experiences a sig- 
nificant loss in angular momentum when it meets off-axis 
shocks during galaxy rotation, this does not necessarily trans- 
late into the mass inflow all the way to the galaxy center. In 
models with c s < 10 km s -1 , a narrow nuclear ring formed 
by gas with X2 -orbits inhibits further inflows of the gas, re- 
sulting in the mass inflow rate M through the inner bound- 
ary less than 0.01 M yr -1 , regardless of the BH mass. The 
mass inflow rates are greatly enhanced to M > 0.01 M yr -1 
in models with M BH < 4 x 10 7 M and c s > 15 km s -1 or 
^bh = 4 x 10 8 M and c s = 20 km s -1 for different reasons 
depending on the BH mass. When M BH < 4 x 10 7 M , the 
gas orbits are affected by thermal pressure and some gas in 
the ring can take highly eccentric orbits, directly falling into 
the inner boundary. In models with M B h = 4 x 10 8 M , on 
the other hand, the gas orbits near the center are more or less 
circular, but the density in the nuclear regions is enhanced 
greatly because of strong nuclear spirals, increasing M. 

6.2. Discussion 

Nuclear rings play an important role in evolution of 
barred galaxies by p roviding sites of active star formation 
near the ce nters ( e. g..lButall l986: Garcia-Barr eto et all 1991: 
iBarth et al.1 119951: iMaoz et al.1 120011: iMazzuca et al.1 l2008h . 
Rings certainly consist of gas that migrates from outer parts 
inward by losing angular momentum at off-axis shocks, but 
what stops further migration to form a ring remains a mat- 
ter of debate. A s a trapping mechanism of the ring material, 
ICombesI £l996) proposed the non-axisymmetric bar toque that 
forces gas to accumulate between two ILRs or at a single ILR 
depending on the shap e of the gr avitational potential (see also 
lButa&Combeslll996l) . while feegan & Teubenl (I2003L 12004ft 
favored the gas transitions from xi- to X2 -orbits rather than or- 
bital resonances. However, our numerical results show that 
a ring is formed at early time because of the centrifugal bar- 
rier that the migrating material feels. Later on, nuclear rings 
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slowly shrink in size as gas with lower angular momentum 
gas is continuously added. Although nuclear rings are located 
in between two ILRs in models with no BH, this is a coinci- 
dence. Models with a central BH show that the specific ring 
positions are insensitive to the shape of the Q-k/2 curves, 
suggesting that nuclear ring formation is not a consequence of 
the orbital resonances. In fact, the gas flows that produce the 
ring are not in force balance and have a large radial velocity, 
so that the concept o f resonances and ILRs is not applicable 
to nuclear rings (e.g jRegan & Te uben 2003). In addition, the 
notion of %i -orbits as the gas trapping locations is meaningful 
only for small c s H\ When the sound speed is large, thermal 
pressure at the contact points between the off-axis shocks and 
nuclear ring causes the gas orbits to deviate from -orbits 
considerably. 

The results of our simulations suggest that not all barred- 
galaxies possess nuclear spirals at their centers. Long-lasing 
nuclear spirals exist only when either the sound speed is 
small or the BH mass is large: they do not survive in 
models with both large c s and small M B h- Two com- 
mon views regarding the nature of nuclear spirals are low- 
ampl itude density waves and strong gaseous shocks (see 
e.g., Enslmaier & Shlosman 2000; Macieiewski et al. 2002; 
Maci ejewskil l2004allbl : lAnn & Thakurl 120051: iThakur et alJ 
2009). And, our simulations indeed show that they are either 
tightly-wound density waves or shocks when the pith angle 
is relatively large (i p > 6°). One may speculate that nuclear 
spirals are more likely to be shocks rather than density waves 
when c s is small. In contrast to this prediction, however, nu- 
clear spirals in models with Mbh = 4 x 10 8 M are density 
waves when c s is small and shocks when c s is large. This is 
of course because as c s increases, waves in nuclear regions 
tend to be more open (with smaller \k\) in the beginning, and 
they are subsequently supplied with more gas from the rings 
as they grow. This i s entirely consistent with the results of 
lAnn & Thakur ( 2005) who used SPH simulations to show that 
nuclear spirals are supported by shocks when c s £ 15 km s -1 
in models with a massive BH. We note however that weak 
trailing spirals seen in our Model csl0bh8 are a bsent in Model 
M2 (c s = 10 km s" 1 and M BH = 4 x 10 8 M ) of lAnn & Thakun 
( 2005), which is presumably due to an insufficient number of 
particles to resolve nuclear spirals in their SPH simulations. 

Of 12 models with differing c s and M B h(0) that we have 
considered, only 1 model possesses leading spirals, suggest- 
ing that galaxies with leading nuclear spirals would be very 
uncommon in nature. To our knowledge, only two galax- 
ies, NGC 1241 and NGC 6902, ar e known to possess lead- 
ing features in the nuclear regions (iDiaz et al.ll2003t lGrosb0ll 
20030. Based on our simulations, the existence of leading 
spirals at centers requires two stringent conditions: (1) the 
gas should be dynamically cold enough to protect nuclear 
spirals from nuclear rings and (2) there should be a wide 
range of radii with d(Q- K,/2)/dR < in the central parts, 
which can be easily accomplished when there are two ILRs 



(or without a strong central mass concentration). The sec- 
ond condition is consistent with the linear theory that pre- 
dicts sho rt leading waves pr opagating outward from the in- 
ner ILR (lMaciejewski l l2004ah. The facts that NGC 6902 is 
a barred- spiral galaxy (lLaurikainen et al.ll2004f> and does not 
show significant X-ray e missions indicative of AGN activities 
(Desroches & Ho 2009) are not inconsistent with the second 
requirement for the existence of leading nuclear spirals. Since 
NGC 1241 is a Seyfert 2 galaxy with an estimated BH mass 
of log(M BH / M ) = 7.46 dBian & Gull2007h . however, the nu- 
clear star-forming regions in this galaxy are unlikely to be 
associated with gaseous nuclear spirals studied in this work. 

Finally, we discuss the mass inflow rates derived in 
our models in regard to powering AGNs in Seyfert galax- 
ies. The mass accretion rate is often measured by 
the Eddington ratio defined by A = Lboi/^Edd = 4.5 x 
10- 2 (e/0.1)(M/10- 2 M yr^XAfBH/lO 7 M ) _1 , where L bo i 
and L E dd denote the bolometric and Eddington luminosities 
of an AGN and e is the mass-to-energy conversion efficiency 
of the accreted material. Observations indicate that A < 0.1 
for cl assical Seyfert 1 galaxies with broad iron emis sion lines 
(e.g., iMeyer-H ofmeister & Meyej] 120111: see also iPetersonl 
11997b . wh ile A ~ 10" 3 for low-luminosity Seyfert 1 AGNs 
(e.g., iHol |2008[) . In our numerical models, the mass in- 
flow rates are larger for models with smaller M rh and larger 
c s . Taking e « 0.1 (e.g.. lYu & Tremainei r2002). the mass in- 
flow rates amount to A <, 10~ 4 for M B h = 4 x 10 8 M and 
c s < 15 km s" 1 , A - 1(T 3 for M BH = 4 x 10 8 M and c s = 
20 km s" 1 or for M B h = 4 x 10 7 M and c s = 10 km s~\ and 
A~0.02-0.4forM BH = 4x 10 7 M andc, > 15kms _1 . For 
classical Seyfert galaxies with strong bars, this suggests that 
the masses of central BHs are likely to be less than 10 8 M , 
which appears to be consistent with the measured values from 
the relation between the BH masses and the velo city dis- 
persions of stellar bulges (e.g., Watabe et al. 2008) and re- 
verberation m apping techniques (e.g., iGultekin e t al. 2009; 
iDenney et al.ll201Ci|) . Of course, this result may depend on 
many factors such as the axis ratio and strength of the bar, 
presence of self-gravity and magnetic fields, gas cooling and 
heating, turbulence, etc., all of which would affect gas dy- 
namic significantly. Extending the present work to include 
these physical ingredients would be an important direction of 
future research. 
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APPENDIX 
THE GALAXY MODEL 

In this Appendix we describe the gravitational potential of the model galaxy that maintains the rotation of the (non- self- 
gravitating) gas disk. The potential is comprised of four component: stellar disk, bulge, bar, and BH. For the disk, we take a 

7 The models considered by| Regan & Teubej GSElEOOl) had the sound ft [s uncertain h()w these ^ features ^ tQ gaseQus 
speed fixed to c s - 5 km s . spirals studied in this paper. 

8 Leading arms in NGC 1241 are detected by Pao; emissions tracing young 
stars, while those in NGC 6902 are observed in the K 7 band tracing old popu- 
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^) = ^rr(l + ^ , (AD 



Kuzmin-Toomre model with surface density 

27tGR V 

where vo and Ro are constants (iKuzminll 1 9561: iToomrdl 1 963h . The corresponding gravitational potential at the disk midplane is 

* disk = --±^, (A2) 

(Binnev & Tremaine 2008). We take vo = 260 km s -1 and Ro = 14.1 kpc in our simulations0 The total disk mass is M d i sk = 
v%Ro/G = 2.2 x lO n M . 
For the bulge, we use a modified Hubble profile with volume density 



R z 



-3/2 



p(«) = Pbui( 1 + ^ ) , (A3) 



and the potential 



* bul = ^ln - + Wl + ^ , (A4) 



where pbui and Rt are the central density and the characteristic size of the bulge, respectively. We take pbui = 2.4 x 10 10 M kpc 3 
and Rb = 0.33 kpc, with th e correspondin g bulge mass of M bu i = 2.8 x 10 10 M within R = 6 kpc. 
The bar is modeled by a Ferrers (1887) ellipsoid with volume density 



-{ 



Pbar(l-gT forg<l, (A5) 

elsewhere, 



where pbar is the bar central density, g 2 = y 2 /a 2 + (x 2 +z 2 )/b 2 with (x,y,z) being the Cartesian coordinates, and a and b are the 
semimajor and semiminor axes of the bar, respectively. The exponent n measures the central concentration of the bar density 
distribution. In this work, we fix the bar parameters to n = 1, a = 5 kpc, b = 2 kpc, and pbar = 4.5 x 10 8 M kpc -3 for all models. 
The total mass of the bar is then M bar = 2 2n+3 7rab 2 p h3iY T(n+ l)r(>T + 2)/r(2ft + 4) = 1.5 x 1O 1O M and the bar quadrupole moment 
is Q m = Mbar<2 2 [l - (b/a) 2 ]/(5 + 2ri) = 4.5 x 10 10 M kpc 2 , corresponding to a strong bar. When n = 1, the bar potential is given 
explicitly as 



$ bar (x, y,z) = ~ ^ Ga f Pbar [ WW +*Vw 20 o + 2v 2 W 110 - 2W U 



2 L - -WWW • " V" --UV, ■ - , •■ 11W -.^100) 

V(/Wo2o + 2z 2 Woii-2W io) 

+z 2 (z 2 W 02 + 2x 2 W 10 i -2Wooi)], (A6) 

where the coefficients WJ-^'s are defined in Pfenniger (1984). As the bar pattern speed, we choose = 33 km s -1 kpc -1 which 
places the corotation resonance at Rqr = 6 kpc. 
Finally, for a central BH with mass M B h, we use a Plummer potential 

*«» = - CAT) 

with the softening radius 7^ = 1 pc. We take M B h = 0, 4 x 10 7 M , and 4 x 10 8 M to study the effects of the central mass 
concentration on the bar substructures. With Mbh <C Mdisk, Mbui, M^, the BH affects the rotation curve only in the very central 
regions (with R < 0.5 kpc). Figures [T] and [2] plot the rotational velocity for models with Mbh = 4 x 10 8 M and the angular 
frequency curves for all models, respectively. 
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